clear; clc; close all;

%% 参数设定
% 初始条件：[x1, y1, x2, y2, x3, y3, phi]
X0 = [-2, 0, -1, -1, -1.12, 0.1, 0.3];

% HR 神经元1 参数
params.a1 = 1.5;    % 激活参数：影响膜电位与恢复变量之间的相互作用
params.b1 = 3;      % 非线性抑制系数：调节膜电位的非线性响应
params.c1 = 1;      % 恢复变量：决定恢复变量的响应特性
params.d1 = 5;      % 对 HR 神经元1的膜电位平方项产生抑制作用
params.i1 = 0.45;   % 外部输入电流/偏置，为 HR 神经元1提供外部激励

% FN 神经元 参数（FitzHugh-Nagumo 模型）
params.a2 = 0.77; 
params.b2 = 0.3;
params.c2 = 0.8;
params.i2 = 0;

% HR 神经元2 参数
params.a3 = 0.9;
params.b3 = 3;
params.c3 = 1;
params.d3 = 5;
params.i3 = 0.5;

% 耦合系数（描述神经元之间的连接强度）
params.m12 = 0.785; % FN - HR1 的耦合强度
params.m21 = 0.6;  % HR1 - FN 神经元的耦合强度
params.m23 = 0.785; % FN - HR2 的耦合强度
params.m32 = 1.01;  % HR2 - FN 神经元的耦合强度

% 分岔参数 ε 的取值范围
eps_range = linspace(0.6, 2.4, 250);

% 积分参数
T_total = 1000;   % 总积分时间
T_trans = 500;    % 舍去瞬态时间

%% 初始化存储数组
eps_cell = cell(length(eps_range), 1);
x1_cell  = cell(length(eps_range), 1);

%% 并行计算每个 ε 值下的局部极值点
parfor idx = 1:length(eps_range)
    eps_val = eps_range(idx);
    
    % 定义 ODE 系统句柄，调用封装函数 hrfnhr
    odefun = @(t, X) hrfnhr(t, X, eps_val, params);
    
    % 使用内置 ode45 求解器积分
    tspan = [0 T_total];
    options = odeset('RelTol',1e-6,'AbsTol',1e-8);
    [t, X_sol] = ode45(odefun, tspan, X0, options);
    
    % 舍去瞬态部分，只保留 t > T_trans 的数据
    ind = t > T_trans;
    t_post = t(ind);
    X_post = X_sol(ind, :);
    x1_signal = X_post(:, 1);
    
    % 提取局部极值点：局部最大值
    [maxVals, ~] = findpeaks(x1_signal, t_post);
    % 局部最小值（对信号取负后找峰值，再取负恢复）
    [minVals, ~] = findpeaks(-x1_signal, t_post);
    minVals = -minVals;
    
    extVals = sort([maxVals; minVals]);
    if isempty(extVals)
        extVals = x1_signal(end-9:end); % 若无局部极值，取后期最后 10 个点
    end
    
    % 存储当前 eps_val 的结果
    eps_cell{idx} = repmat(eps_val, length(extVals), 1);
    x1_cell{idx}  = extVals;
end

%% 将 cell 数组合并成向量
eps_plot = vertcat(eps_cell{:});
x1_plot  = vertcat(x1_cell{:});

%% 绘制分岔图
figure;
scatter(eps_plot, x1_plot, 1.5, 'k', 'filled');
xlabel('\epsilon', 'FontSize', 12);
ylabel('x_1', 'FontSize', 12);
title('x_1 关于 \epsilon分岔图', 'FontSize', 14);
ylim([0.7 1.1]); 
set(gca, 'FontSize', 12);